#ifndef _FDM_H
#define _FDM_H

/**
 * @file FDM.h
 * @author your name (you@domain.com)
 * @brief 有限差分
 * @version 0.1
 * @date 2022-05-31
 * 
 * @copyright Copyright (c) 2022
 * 
 */

#include <iostream>
#include <vector>
#include <functional>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Sparse>
#include <cmath>
#include <set>

using std::vector;
using std::function;
using std::set;
using namespace Eigen;
using Tri = Eigen::Triplet<double>;
using SpMat = Eigen::SparseMatrix<double>;

class FDM_DIM1
{
private:
    size_t nx;      //在x方向上的分段数
    double x_lhs;   //左端点的坐标
    double x_rhs;   //右端点的坐标
    function<double(double)> f_rhs;      //possion方程右端项
    function<double(double)> u_realvalue;   //方程解析解，即边界解
    SpMat A_lhs;    //拼装出的矩阵
    VectorXd F_rhs;     //拼装出方程的右端项
    set<size_t> boundary;       //边界点
public:
    FDM_DIM1() = default;
    FDM_DIM1(size_t _nx, double _x_lhs, double _x_rhs, function<double(double)> _f_rhs, function<double(double)> _u_realvalue):
            nx(_nx), x_lhs(_x_lhs),x_rhs(_x_rhs),f_rhs(_f_rhs), u_realvalue(_u_realvalue) { };
    SpMat& get_matrix() { return A_lhs; };
    VectorXd &get_Frhs() { return F_rhs; };
    set<size_t> &get_boundary() { return boundary;};
    void build();
};

void FDM_DIM1::build()
{
    double h = 1.0 *(x_rhs - x_lhs) / nx ;
    size_t n_vec = nx + 1;
    VectorXd F = VectorXd::Zero(n_vec);
    size_t i = 0 ;
    F[0] = u_realvalue(x_lhs) /(h*h);
    vector<Tri> Trilist;
    Trilist.push_back(Tri(0,0,1.0));
    Trilist.push_back(Tri(0,1,0));
    i++;
    for( ; i < n_vec - 1; i++)
    {
        Trilist.push_back(Tri(i,i-1,-1.0));
        Trilist.push_back(Tri(i,i,2.0));
        Trilist.push_back(Tri(i,i+1,-1.0));
        F[i] = f_rhs(x_lhs + h * i);
    }
    Trilist.push_back(Tri(n_vec-1, n_vec -1, 1.0));
    Trilist.push_back(Tri(n_vec-1, n_vec -2, 0.0));
    F[n_vec - 1] = u_realvalue(x_rhs) /(h*h);
    SpMat A(n_vec, n_vec);
    A.setFromTriplets(Trilist.begin(), Trilist.end());
    A.makeCompressed();
    this->A_lhs = A /(h*h);
    this-> F_rhs = F ;

    set<size_t> bound;
    bound.insert(0);
    bound.insert(n_vec - 1);
    this->boundary = bound ;
};


#endif